Heterogeneous slow dynamics in a two dimensional doped classical antiferromagnet 
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We introduce a lattice model for a classical doped two dimensional antiferromagnet which has 
no quenched disorder, yet displays slow dynamics similar to those observed in supercooled liquids. 
We calculate two-time spatial and spin correlations via Monte Carlo simulations and find that for 
sufficiently low temperatures, there is anomalous diffusion and stretched-exponential relaxation of 
spin correlations. The relaxation times associated with spin correlations and diffusion both diverge 
at low temperatures in a sub-Arrhenius fashion if the fit is done over a large temperature-window 
or an Arrhenius fashion if only low temperatures are considered. We find evidence of spatially 
heterogeneous dynamics, in which vacancies created by changes in occupation facilitate spin flips 
on neighbouring sites. We find violations of the Stokes-Einstein relation and Debye-Stokes-Einstein 
relation and show that the probability distributions of local spatial correlations indicate fast and 
slow populations of sites, and local spin correlations indicate a wide distribution of relaxation times, 
similar to observations in other glassy systems with and without quenched disorder. 

PACS numbers: 05.20.-y, 75.10.Hk, 75.10.Nr 



I. INTRODUCTION 

Doped two-dimensional antifcrromagnets have at- 
tracted considerable attention in the past two decades, 
due to their relevance as a model to describe high tem- 
perature superconductivity in the cupratesi^ In addition 
to superconductivity, there have been many other phases 
observed in the cuprates, and there has been particular 
recent interest in the spin glass observed at low temper- 
atures and intermediate doping) 3 ^ 4 * 5 ^ and the checker- 
board ordered electronic state in the high-T c cuprate 
P^S^CaCi^Os+d" (Bi2212) and in hole-doped copper 
oxides. 7,8 Whilst quantum effects are very important in 
the cuprates, it is also of interest to study whether glassy 
phases and slow dynamics are present in simple classical 
doped two dimensional antiferromagnets, and this is the 
subject we address here. 

In this paper we construct the simplest such model 
that we can find that displays glassy dynamics - it has 
nearest neighbour interactions only and the spins are on 
a square lattice, so there is no geometric frustration. The 
only source of frustration is that there are mobile holes, 
and these mean that local relaxation to Neel order leads 
to local frustration. 

Whilst the model we construct is principally that of 
an antiferromagnet, it falls within the general class of 
models which show slow dynamics without quenched dis- 
order. There has been intense activity in the past few 
years investigating models with this behaviour, espe- 
cially kinetically constrained models, which have triv- 
ial Hamiltonians but more complicated rules governing 
their dynamics, 9 ! 10 ! 11 : 12 ! 13 ! 14 ! 15 : 16 ! 17 ! 18 These models have 
been investigated with the aim of shedding light on the 
local structure of slow dynamics in glasses, dynamical 
heterogeneitiesii2iSS*2i The approach here is complemen- 



tary to the work on kinetically constrained models, as 
we study a model which has a more complicated Hamil- 
tonian, but relatively simple dynamics. The model we 
study also has some similarities to frustrated lattice 



,22,23 



and the hard square lattice gas. 



gasc: 

We find that this model displays slow dynamics at low 
temperatures, both in the diffusion of particles and in 
the relaxation of spins. For the values of the param- 
eters used (temperature, doping and ratio between the 
antiferromagnetic and nearest-neighbour repulsion) and 
the linear sizes analyzed the system reaches a station- 
ary state after a transient. This state is characterized by 
antiferromagnetic or checkerboard order. In the station- 
ary state we find unusual temperature dependence of the 
relaxation times for both diffusion and spin relaxation, 
that can be fit by a sub-Arrhenius temperature depen- 
dence over the entire temperature range that we study 
(although there is Arrhenius temperature dependence at 
low temperatures). We also find evidence of spatially 
heterogeneous dynamics that lead to a breakdown of the 
relation between diffusive and rotational motion dictated 
by the Debye-Stokes-Einstein law (where we consider the 
spin degree of freedom to model rotations). The distri- 
butions of local correlations measured at different time 
differences are stationary, and hence trivially scale with 
the value of the global correlation^ The form of these 
distributions is interesting in that they suggest that for 
spatial correlations there are two populations of sites, one 
with fast, and the other with slow dynamics, and a wide 
range of spin relaxation times. 

The paper is structured as follows. In Sec. |n] we in- 
troduce the model that we study and the one-time and 
two-time quantities that we calculate with Monte Carlo 
simulations. In Sec. II I II we display our results for the 
phase diagram and the correlation functions that we de- 
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fined in Sec. [n] We show evidence of time-scales that 
diverge at low temperatures, and find that both the dif- 
fusive motion and spin flip dynamics are spatially het- 
erogeneous. We also investigate the distributions of local 
two-time quantities. In Sec. IIVI we discuss how our re- 
sults indicate a breakdown in the Debye-Stokes-Einstein 
relation. Finally, in Sec.^we summarize and discuss our 
results. 



neighbouring sites or to flip its spini^ We only allow at- 
tempts to move a particle to an unoccupied site (with 
equal probabilities for the unoccupied sites) iSL The at- 
tempt is then accepted with Boltzmann probability. In 
each of the particle move and spin flip steps, detailed bal- 
ance is respected. The Monte Carlo time-unit is equal to 
the attempt of N updates (motion of particles or spin 
flips). 



II. MODEL 



Motion of holes 



The Hamiltonian for the model is 



n 



riiUj (V + JSiSj) , 



(1) 



where ni = 0, 1 is a density variable indicating occupa- 
tion of a site, and Si = ±1 is an Ising spin attached to 
each particle (i.e. when a particle moves site, so does 
its spin). There is a nearest neighbour repulsion with 
magnitude V, and a nearest neighbour antifcrromagnetic 
interaction with magnitude J. The particles have a hard- 
core constraint so that there is no double occupancy of 
sites. We are interested in the limit where J/V < 1 (in 
the limit where J/V > 1, we find phase segregation of an- 
tiferromagnetic domains and regions with no particles). 
We study this model in the canonical ensemble by fix- 
ing the number of particles on a two dimensional square 
lattice and define the hole concentration x = 1 — N/L 2 , 
where N is the number of particles and L is the lattice 
size. This is essentially the classical, Ising limit of the 
t — J — V model studied by Kivelson and Emery and 
co-workers 

We note that the model has some similarities to others 
that have been introduced in the literature, in particular 
the frustrated lattice gas^^ although in that case, the 
spin interactions Jy are randomly drawn from a distribu- 
tion, rather than being of constant sign and magnitude. 
In the limit that J — > oo, this model has some resem- 
blance to the hard square lattice gas24 (although the spin 
degree of freedom means that the dynamics here have a 
different flavour to that model). 

When there is no doping, the ground state is an antifer- 
romagnet with Neel order. This implies that the system 
can be divided into two sublattices, each with ferromag- 
netic order, shifted by one lattice spacing in the x and y 
directions with respect to each other. In this limit, the 
model can be made equivalent to a ferromagnet by apply- 
ing a spin-flip transformation on one sublattice. However, 
once there is doping and holes are free to move around, 
the transformation is no longer applicable and the anti- 
ferromagnet and ferromagnet are distinct. 

We simulate this model using classical Monte Carlo 
simulations and Metropolis dynamics. The model has 
no intrinsic dynamics, hence we impose the following dy- 
namics. At each step we choose with 50 % probabil- 
ity either to attempt to move a particle to one of its 



Whilst naively it might be thought that an individual 
hole would be localized due to leaving a string of over- 
turned spins behind it, holes are able to move diagonally 
in an antifcrromagnetic background at no energy cost by 
hopping around one and a half loops of a plaquette. 30 
The energy of the final state is the same as that of the 
initial state, however this is an activated process with 
activation energy 12 Jj2£ as the hole flips spins into an 
unfavourable energy configuration on its first lap of the 
plaquette then lowers the energy to that of the original 
state as it goes round for another half a lap. Similarly, a 
pair of holes that are on opposite corners of a plaquette 
can move along the diagonal that they lie on, although 
the activation for this process can involve V as well as J. 
At low doping, where holes interact very rarely, this im- 
plies that there are a lot of states with the same energy, 
yet with large barriers between them, which naturally 
leads to slow dynamics. 



A. Quantities Calculated 

We calculate quantities that should demonstrate spin 
and hole ordering and illustrate the dynamics of the spins 
and holes. The one-time quantities that indicate order- 
ing are the staggered magnetization that indicates Neel 
order of the spins and a checkerboard order parameter 
that indicates ordering of the holes. As has been dis- 
cussed extensively in the literature (see e.g. Refs. I3lll3^1 . 
glassy dynamics are often seen more easily in two-time 
quantities if one-time quantities reach a limit. In most 
of the cases we consider the one-time quantities saturate 
for high temperatures, but can saturate quite slowly as 
the temperature is decreased. We calculate the two-time 
correlations for spins and mean square displacement for 
both high and low temperatures. 



1. One time quantities 

We calculate two one-time quantities. These are the 
staggered magnetization 



M„ = 



1 E(- 



.r,i/ 



x.y 



(2) 
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and a checkerboard order parameter 



B. Parameters 



Mr 



N ^ 



(3) 



where x and y are the co-ordinates of a given lattice site 
and S Xi y and n x ,y are the spin and occupation of the site 
i — (x,y) respectively. We calculate these quantities as a 
function of temperature to determine the phase diagram, 
which is shown in Fig. 



In all our runs we used random initial conditions. Most 
of the data presented are for L = 30 though we also used 
larger systems L = 40, and 50, to determine the phase 
diagram. We set J/V = 0.2 unless otherwise specified. 
In our calculations we mainly use a waiting time of t w — 
3200 Monte Carlo steps (MCs) (we found that we had 
essentially identical results for a shorter waiting time of 
t w = 320 MCs, and a longer waiting time of t w = 32000 
MCs). We accumulated data for up to 10 8 MCs after the 
waiting time. 



2. Two time quantities 

We calculate several two-time correlation functions, for 
both spatial and spin correlations. These are: the mean 
square displacement 

^M-H^EIM*)- 1- "^)) 2 !' ( 4 ) 

a 

where r a is the position of the a th particle; a spin auto- 
correlation function 

C(t, t w ) = —^2s a (t)s a (t w ), (5) 

a 

where s a is the spin of the a th particle; and a constant- 
site spin correlation function 

Clocal{t,t w ) = Si(t)Si(t w ), (6) 

i 

where Si is the spin at site i, and we note that Si 
may not be the same spin at times t and t w (although 
at low enough temperatures they coincide unless t and 
t w are very widely spaced). We expect the long time 
limits of C(t,t w ) and Ci oca i(t,t w ) to have different be- 
haviour if there is Neel order present in the system. As 
t — > oo, we expect C(t,t w ) — > 0, whereas Ci oca i(t,t w ) — > 
M s (t)M s (t w )~ The arguments for each of these lim- 
its are as follows. If there is Neel order, then the lat- 
tice can be divided into two sublattices, each of which 
has an opposite spin orientation. Consider the spins on 
one sublattice, at long times t, diffusion will mean that 
half are on the other sublattice, and will have their spin 
flipped relative to its orientation at time t w . The con- 
tributions of the set of spins on the two sublattices to 
C(t, t w ) will have opposite signs and equal magnitude and 
hence we expect C(t, t w ) — > 0. In contrast, if we consider 
the same-site correlation, then if there is no flipping of 
the antiferromagnetic order, Ci oca i(t, t w ) should decay to 
M s (t)M s {t w ) reflecting the average value of the staggered 
magnetization on the site at the two times. We mainly fo- 
cus on C(t, t w ) but show that Ci oca i(t, t w )- M s (t)M s (t w ) 
has qualitatively similar time dependence to C(t,t w ). 



III. RESULTS 

A. One-time quantities 

We first display results for one-time quantities which 
demonstrate that these saturate relatively quickly at high 
temperatures, but the saturation can be quite slow as the 
temperature is lowered at low values of x. The staggered 
magnetization at five different temperatures is shown as 
a function of time in Fig. ^for x = 0.1 and J/V = 0.2, 
averaged over 50 samples of size L — 30. For this choice 
of parameters, the configurations have local Neel order- 
ing at relatively short times, as shown in Fig. [21 for one 
sample at low temperature with one large domain. How- 
ever, there can be domain walls that move slowly, and 
lead to a relatively slow saturation of the staggered mag- 
netization. The two-time correlations defined in Eq. (jSJ) 
do not appear to be particularly sensitive to the presence 
of antiferromagnetic domain walls. 
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FIG. 1: (Color online) Staggered magnetization as a function 
of time averaged over 50 samples with L = 30 for x = 0.1, 
J/V = 0.2 and temperatures T/J = 1.2, 1.0, 0.8, 0.6, and 
0.4. 

At larger values of x, there is checkerboard order in 
addition to Neel order, and in Fig. |3 we show the time- 
evolution of the checkerboard order parameter M c and 
staggered magnetization M s averaged over many sam- 
ples for L = 30 with x — 0.35 at four different temper- 
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FIG. 2: Illustration of the state of the system after 400 MCs in 
a sample with L = 30 for x = 0.1, J/V = 0.2, and T/J = 0.1. 
Up spins are indicated by black squares and down spins are 
indicated by open squares. 
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atures. Figure 0] shows the configuration reached by a 
low temperature sample at the same doping level after 
25600 MCs. The results in the two figures appear to 
be consistent with the coexistence of antiferromagnetic 
and checkerboard order at low temperatures (the non- 
zero value of M s at T > 0.8 in Fig. b) is a finite size 
effect), and Fig. 0] appears consistent with the possibil- 
ity of phase separation. We note that the time-scale for 
the onset of checkerboard order grows as temperature is 
decreased so that for T = 0.4, the checkerboard order 
has not saturated in the time-window considered, even 
though at higher temperatures in Fig. [3] a) the value of 
the order parameter saturates at around M c ~ 0.47. 

We show the phase diagram in Fig. [S] and indicate the 
regions in which there is Neel ordering and checkerboard 
ordering. The boundary of the region of Neel order, the 
Neel temperature, Tjy, is found by calculating M a at 
L = 30, 40, and 50 and extrapolating to where M s van- 
ishes. We also find that for x = 0, T N = 2.29 ± 0.05 J, 
in agreement with the exact value of Tn — 2.27 J for the 
two-dimensional Ising antiferromagnet. We use the same 
procedure with M c to calculate the checkerboard order- 
ing temperature T c (,. The third temperature shown on 
the phase diagram is T ne , which is the temperature below 
which we see stretched exponential spin relaxations, and 
is discussed further in Sec. IIII Bl From Fig. it is clear 
that Neel order disappears as x is increased towards 0.5, 
and checkerboard order grows for x > 0.2. For immobile 
holes, one expects the antiferromagnetism to disappear 
at the percolation threshold. The percolation threshold 



FIG. 3: (Color online) a) Checkerboard order parameter and 
b) staggered magnetization as a function of time in Monte 
Carlo steps averaged at least 100 samples with L — 30 for 
x = 0.35, T/J = 0.4, 0.6, 0.8, and 1.0, and J/V = 0.2. 



for bond dilution in two dimensions is x c = 0.5*2* and for 
site dilution it is x c = OAlm^^ For any one snapshot of 
the spin-hole configurations, there is no percolating clus- 
ter for x > 0.41 in the thermodynamic limit, and hence 
we would expect that x c = 0.41. Whilst we have not 
investigated this question in great detail, our results are 
consistent with antiferromagnetism vanishing at the site 
dilution threshold. The slow dynamics typically occur at 
temperatures lower than the Neel temperature, although 
we note that T ne > Tjv for x = 0.4, in the region where 
T c b > T ne > T^. 

Changing the ratio of J/V whilst keeping it less than 
1 appears to have little effect on the Neel order, but the 
checkerboard transition temperature decreases rapidly as 
J/V increases. For instance, if J/V — 0.5 then we find 
that T c b/J is reduced to 1.3 at x — 0.5, compared to 
T cb /J = 3.8 for J/V = 0.2, as shown in Fig. This 
appears to indicate that T c b ~ V at low J/V, whilst 
Tn ~ J; we have not explored the dependence of T ne on 
J and V. 



5 



up spins 
down spins 



■□ □ 



■ ■ ■□ □■ 

■ ■■□□[ 

■ ■ ■□ □ i 



□ ± 
□ □ 



□ ■□ 

1 □ c 

□ □ 



□■□ □ ■ c 
□ □ ■ □■ 



□ ■□ ■ 

■ □ □ 

nan 

■ □ □ 



□ □ 
□ □ ■ 



□ □■ 

■ □ □■ 

□ ■□ c 



□ □ 
■ □ c 

□ □■ 



30 

25 

20 

15 

10 

5 





5 10 15 20 25 30 



FIG. 4: A configuration with checkerboard hole-ordering, and 
Neel ordering for x = 0.35, L = 30, J/V = 0.2, and T/J = 0.1 
after 25600 MCs. 




FIG. 5: (Color online) Phase diagram, showing Neel temper- 
ature Tjv, checkerboard ordering temperature T c t, and onset 
temperature of non-exponential spin relaxation, T„ e as a func- 
tion of x and T/J, for J/V = 0.2. The antiferromagnetic (AF) 
and checkerboard (CB) regions are marked, as is the region of 
apparent coexistence of antiferromagnetic and checkerboard 
order (AF + CB). 



B. Two-time quantities 

We show the spatial and spin correlation functions as a 
function of t — t w at several different temperatures below. 
All data shown here for two-time correlations was taken 
in L = 30 systems - the data shown is only for one ther- 
mal history, although it was checked that changing ther- 
mal history does not quantitatively change the results or 
fits. We note that there is some waiting time dependence 
of the results at low temperatures, although this is gen- 
erally only for very short waiting times (t w < 320 MCs). 



Hence for the data we show here, we actually find that 
the two-time quantities depend only on time differences 
for the waiting times we use, so D(t, t w ) ~ D(t — t w ) and 
C(t, t w ) — C(t tw)- 

Our analysis of the data follows the following scheme. 
We attempt to fit 



D(t,t w ) = ((t-t w )/r r ) a , 



(7) 



and 



C(t,t w ) = A.e-^'^/^, (8) 
Ci ocal {t,t w ) = Aie-W-W'V + B, (9) 

where r r is the relaxation time for diffusion and t s is the 
relaxation time for spin correlations. Apart from very 
short times and very long times (where the system size 
becomes important for D, and spin correlations die out 
for C), these forms work very well (note that A\ — > 1, 
A2 — > 0, and B — > 1 at low temperatures and low dop- 
ing). We use the fits to the two-time spin correlation to 
fit T ne , the highest temperature at which one observes 
non-exponential relaxation, i.e. fi < 1. This tempera- 
ture scale is not associated with a change in any order 
parameter that we measure, but allows us to determine 
the boundary of the region in the phase diagram where 
slow dynamics are observed. However, we do note that 
the equilibration of one-time quantities at T < T ne is 
generally much slower than for T > T ne . 

In the figures below, we only show data from x = 0.1 
unless indicated, but it is representative of the behaviour 
seen at other levels of doping, where we observe similar 
qualitative behaviour, with only quantitative differences. 
In Fig. we show the mean square displacement as a 
function of time for x — 0.1 for J/V = 0.2. 
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FIG. 6: (Color online) Spatial correlations for x = 0.1, L = 
30, and J/V = 0.2 for temperatures T/J = 0.4, 0.8, 1.2 
and 1.8. Note that at the lowest temperature the slope of 
the curve changes, indicating anomalous diffusion. We also 
indicate l\ = l/4s = 2.5 and the fits to D(t, t w ) at T = 0.4 for 
D(t,t w ) < i| and D(t,t w ) > 1%. Note also that the saturation 
at long times and high temperatures is a finite size effect. 
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Note that at the lowest temperature in Fig. there is 
a signature of anomalous diffusion (i.e. a < 1). This 
is present for mean square displacements that are less 
than l 2 x — l/ix, which is the square of the lengthscale 
equal to half the average distance between holes. We 
show a more striking example of anomalous diffusion in 
Fig. [7| for x = 0.02, where the lengthscale l 2 x is much 
larger than in Fig. In general we only see anomalous 
diffusion at temperatures T < T ne and D(t,t w ) < I 2 ., 
however the temperature where we first observe anoma- 
lous diffusion does not appear to be related to T ne , which 
is when the spin auto-correlations show non-exponential 
relaxation. We note that anomalous diffusion has been 
observed in several experimental glassy systemsj^LiS and 
can be understood as "caged" motion on short length- 
scales and timescales, with usual diffusion taking over at 
longer timescales and lengthscales when particles have es- 
caped from their cages. Here, the lengthscale for caging 
is set by the mean distance between vacancies, l x . 
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FIG. 7: Spatial correlations for x = 0.02, L = 30, and 
J/V = 0.2 for temperatures T/J — 0.4. Anomalous diffu- 
sion is clearly evident. We also indicate 1% = 1/Ax = 12.5 and 
a fit to D(t,t w ) of the form ((t - t w )/r) - 8 . 

In Fig. |H1 we show the spin correlations for x = 0.1 
and J/V — 0.2 for several different temperatures. We 
note that the timescale for C(t,t w ) — > is of the same 
order as the timescale for D(t,t w ) ~ 1, corresponding 
to the average particle having moved one lattice spacing. 
The argument explaining why C(t,t w ) — > as t — > oo 
in Sec. Ill A 21 gives a natural explanation of why these 
timescales should be of the same order, and indicates the 
importance of the interplay between spin and diffusion 
dynamics, as a particle will generally be forced to flip its 
spin as it diffuses. For T/J < 1.0 we show both expo- 
nential and stretched exponential fits to C(t,t w ). The 
stretched exponentials are clearly better fits at low tem- 
peratures. We also performed calculations of Ci oca i(t, t w ) 
for the same parameters as in Fig. [5J and these cor- 
relations are shown in Fig. Whilst Ci oca i(t,t w ) de- 
cays to a different limit than C(t,t w ) it also displays 
stretched-exponential spin relaxations at low tempera- 



tures, although at a slightly lower temperature than the 
T„ e defined by C(t,t w ). 
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FIG. 8: (Color online) Spin auto-correlations for x = 0.1, 
L = 30, and J/V = 0.2, with T/J = 1.2, 1.0, 0.8, 0.6, and 
0.4. For T/J < 1.0 we display both exponential and stretched 
exponential fits to the correlations - in each of these cases, 
the stretched exponential is the better fit. 
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FIG. 9: (Color online) Spin same-site correlations for x = 0.1, 
L = 30, and J/V = 0.2, with T/J = 1.2, 1.0, 0.8, 0.6, and 
0.4. For T/J < 0.6 we display both exponential and stretched 
exponential fits to the correlations - in each of these cases, 
the stretched exponential is the better fit. 

In Fig. ^3 we show the stretched exponential (3 pa- 
rameter as a function of temperature for the spin auto- 
correlations shown in Fig. which clearly indicates that 
the spin correlations do not decay exponentially in time 
at low temperatures. We also performed simulations for 
x = 0, in which we found no decay of the spin correla- 
tions, indicating that the timescales observed here arise 
solely due to doping. 

We extract values of r r and t s from fitting D(t,t w ) 
and C(t,t w ) as described above, and use these to define 
relaxation time-scales that we then fit as a function of 
temperature: 



lr,s exp 



(10) 
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FIG. 10: Stretched exponential /3 parameter as a function of 
temperature for x = 0.1, L = 30, and J/V = 0.2. 



E 



o 



a r 



0.02 10.3 0.58 10.9 0.62 

0.1 8.6 0.63 8.4 0.68 

0.2 7.5 0.68 7.7 0.71 

0.3 6.9 0.76 7.1 0.77 

0.4 6.4 0.77 6.7 0.82 



TABLE I: Eq and a data from fits to r r and r s made using 
Eq. JTHJ 



where "f r ,s, r o' S i Eq' s and fitting parameters. We 

also tried fitting our data to a Vogel-Fulcher form, but 
this did not lead to good fits. We show fits of the form in 
Eq. (|10|) for r r and t s determined for x = 0.1, J/V = 0.2 
in Figs. II II and 1121 In Table [I] we show values of Eq and 
a that we have extracted for r r and r s in fits of the form 
of Eq. Ijl Op for different x. Those relating to r r are Eq 
and a r and those relating to t s are Eq and a s . 

The slowing down with temperature is similar to that 
seen in glass formers. However, there is a difference - 
in glass formers, the exponent a is equal to 1 for strong 
glass formers (Arrhenius) and 2 for fragile glass formers 
(super- Arrhenius)^ whereas the relaxation times found 
here consistently have a < 1, i.e. sub- Arrhenius be- 
haviour if one tries to fit over the entire temperature 
range in which there is a 6 order of magnitude change 
in the relaxation time. We did find however, that if we 
restricted the temperature range, then there is Arrhe- 
nius temperature dependence of the relaxation times at 
low temperatures, as is also indicated in Figs.llllandll2l 
This suggests this model has strong-glass like behaviour 
at low temperatures. 

The diffusion relaxation time, r r , and spin relaxation 
time t s , display purely Arrhenius behaviour at low tem- 
peratures (T < T ne ), and it is possible to get a good fit 
to both r r and r s with an Arrhenius form at high tem- 
peratures T > T ne , but neither fit is good over the entire 
temperature range, as is visible in Figs. Illl andll2l 
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FIG. 11: (Color online) Fit to diffusion relaxation times de- 
termined from spatial correlation functions as a function of 
temperature at x = 0.1 and J/V = 0.2. 
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FIG. 12: (Color online) Fits to spin relaxation times deter- 
mined from spin correlation functions as a function of tem- 
perature at x = 0.1, L = 30 and J/V = 0.2. 



C. Spin configurations 

In Figs. and ^] we show how changes in the con- 
figurations of spins and occupancy occur in the region 
where slow dynamics dominate. The doping levels used 
are x = 0.1 (Fig.HSJ) and x = 0.3 (Fig-EJl, and both fig- 
ures are for low temperatures {T/J = 0.1). It is obvious 
from these plots that the dynamics is quite spatially het- 
erogeneous at this temperature. Spin flips occur in the 
vicinity of holes and are facilitated by changes in occu- 
pation on adjacent sites. There are also regions in which 
there are no changes in spin or occupation during the 
time-window employed. This implies that the dynamics 
are a combination of vacancy motion and spin flips rather 
than being associated with domain wall motion. Note 
that the figures show where there has been a net change 
in occupation or a net spin flip, rather than whether this 
is the only change that has occurred between t and t w . 
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FIG. 13: Illustration of the sites at which there has been a 
change in occupancy or a spin flip between times 400 MCs 
and 25600 MCs for x = 0.1, T/J = 0.1, J/V = 0.2, and 
L = 30. 



D. Distributions of correlations 

In this section we show distributions of local spin cor- 
relations and displacements^ To obtain these distribu- 
tions, we coarse-grain the correlations that give D(t,t w ) 
and C{t,t w ) when averaged over all sites, over a small 
region about each site. More precisely, we define the 
coarse-grained local square displacement D^ s (t,t w ) and 
coarse-grained local spin correlation C!/ S (t,t w ) as 
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D?{t,t w ) = -j (r a (t)-r a (t w )) 2 , (11) 

a£Ai(t) 

C- S (M™) = j s a (t)s a (t w ), (12) 

aeAt(t) 

where A is the area of the box Ai(i), which is the box 
centered on site i, containing the set of particles {a} at 
time t. In all of the histograms shown below, we choose 
A = 3 2 = 9. In Fig.^|we show the distribution of coarse- 
grained local square displacements for various values of 
t — t w . We find that for t and t w with the same global 
mean square displacement, i.e. D{t, t w ), the distributions 
collapse on each other. In Fig. ^3 we plot similar dis- 
tributions for coarse-grained local spin auto-correlations 
and also see a clear scaling when the global correlation 
is the same. We note that for time separations such that 
D(t,t w ) is the same, C(t,t w ) also coincides. This scaling 
of distributions of local quantities with the value of the 
global correlation at a given time has been seen previ- 
ously in investigations of other models with slow dynam- 
ics and no quenched disorder^ 



FIG. 14: The two figures illustrate the sites at which there 
has been a change in occupancy or a spin flip between times 
400 MCs and 3200 MCs for a) and between 400 MCs and 
25600 MCs for b). The parameters are x = 0.3, T/J = 0.1, 
J/V = 0.2, and L = 30. 

The distributions shown in Figs. IT51 and 1161 appear to 
be stationary (i.e. depend only on t — t w rather than 
t and t w separately). Since the distribution is station- 
ary, it implies that its first moment is also stationary, 
and hence the distribution itself scales with the global 
correlation as follows from the arguments below. If the 
correlation and mean squared displacement are station- 
ary and monotonic, then they can be written as (using 
the mean-squared displacement as an example) 

D(t,t w )=d(t-t w ), (13) 
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FIG. 15: (Color online) Distribution of coarse-grained local 
square displacements for various values of t — t w , with x = 0.1, 
T/J = 0.1, J/V = 0.2, averaged over 20 samples of size 
L = 50. We show three pairs of data, which are offset from 
each other along the x-axis, and we show the average value of 
the global mean square displacement D av for each pair. 
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FIG. 16: (Color online) Distribution of coarse-grained local 
spin auto-correlations for various values of t — tw with x — 0.1, 
T/J = 0.1, J/V = 0.2, averaged over 20 samples of size L — 
50. We show the average value of the global spin correlation, 
C av for each pair of samples. 



and for a given value of D one finds the associated t — t y 



sional model with quenched disorder— and in a numerical 
simulation of colloidal gellation. — The local spin correla- 
tions do not have a clear separation into two populations 
of spins, but there is clearly a wide distribution of lo- 
cal relaxation times associated with the large spread in 
correlations. 



IV. STOKES-EINSTEIN RELATION 

The model studied here has several features that are 
reminiscent of a glass-forming liquid. At high temper- 
atures, the mean-square displacements are linear in t, 
showing diffusive behaviour, and the relaxation times 
associated with this diffusion grow very quickly with 
decreasing temperature. These diffusive properties are 
what one would expect for a liquid of hard-core parti- 
cles, and what we study is essentially a lattice model 
of this type of system. If the hard-core particles were 
also allowed to have a rotational degree of freedom, then 
again, the spin of the particles may be regarded as mod- 
elling this physics. These similarities between the model 
we study, and a lattice model for a liquid of hard-core 
particles, inspires us to ask whether well-known proper- 
ties of liquids are shared by the model. In particular, it 
is interesting to ask whether relationships similar to the 
Stokes-Einstein (SE) and Debye-Stokes-Einstein (DSE) 
relations hold for translational and spin correlations re- 
spectively. 

The SE and DSE relations describe the translation and 
rotational motion of a large spherical particle of radius R 
in a hydrodynamic continuum in equilibrium at temper- 
ature T with viscosity rj. The SE relation predicts the 
dependence of the translational diffusion coefficient, T>, 
on T, R and n: 



V = 



QuRrj 



(15) 



T> is defined from the long-time limit of the mean-square 
displacement, ((x(t) ~ (x{t'))) 2 ) = 2dV(t - t'), d is the 
space dimension, and k B is Boltzmann constant. Sim- 
ilarly, the DSE relation predicts the dependence of the 
rotational correlation time, i rot on T, R and r\: 



D~ X \D\ =t-t n 



(14) 



We can use similar arguments for C(t,t w ). 

However, whilst the scaling of the distribution with the 
global correlation is trivial, the form of the distribution is 
not, and indicates spatially heterogeneous dynamics. In 
the distribution of local displacements, it is clear that as 
t — t w grows there are two populations of sites, one with 
fast dynamics, and another with slow dynamics. We note 
that the distributions displaying two populations of sites 
have D(t,t w ) > l 2 ,, implying that one population con- 
sists of particles still trapped in cages of size ~ l x and 
the other is those that have escaped these cages. This is 
reminiscent of behaviour seen recently in a two dimen- 



t. 



4iri]R 3 
3k B T ' 



(16) 



with t rot extracted from the decay of, say, —(s(t) ■ s(t')), 
where n is the dimension of the orientation degree of free- 
dom s. Equations (|15f) and Ijltjfl imply that the product 
T>t m t should not depend on T and r\. Even though these 
relations are derived for a spherical tracer, in normal liq- 
uids they are often satisfied for the translational and ro- 
tational motion of generic probes and the constituents 
themselves within a factor of 2. 

In a supercooled fragile liquid the situation changes. 
While the rotational motion is consistent with Eq. |JT§}, 
the diffusion of small probe molecules j^iiii as well as the 
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self diffusion (■aaSi is much faster than would be expected 
from the viscosity dictated by Eq. i|15|). The T depen- 
dence of T> is not given by T/rj but there is a "transla- 
tional enhancement" meaning that, on average, probes 
translate further and further per rotational correlation 
time as T g is approached from above. For example, rota- 
tional motion in OTP is in agreement with Eq. i|16|l over 
12 decades in viscosity while the deviation in the trans- 
lational diffusion of small probe molecules of a variety of 
shapes, and the OTP molecules themselves, is such that 
the product T>t Iot increases by up to four orders of mag- 
nitude close to T g . These results imply that molecules 
translate without rotating much more than DSE-SE re- 
lations would require. 

The mismatch between the translational and rota- 
tional motion has been observed numerically2ML2£ and 
it can be described assuming the existence of dynamic low 
viscosity regions, 43,49,50 a free-energy landscape based 
model^ within the random first order transition sce- 
nario^ and kinetically facilitated spin modelsiSi^ 

Our results in Sec. lIII 51 where we demonstrate anoma- 
lous diffusion in Figs. [5] and indicate that the SE rela- 
tion breaks down at low temperatures in this model. To 
test whether an analogy to the DSE holds in this model, 
in Fig. El we plot the temperature dependence of the 
product of the translational diffusion coefficient and the 
characteristic time for relaxation of the spin correlation 
for five different values of x. We find that that there 
is strong x dependence of the ratio t s /t t ~ T>t rot at low 
temperatures: for small x, where there is no checkerboard 
ordering, r s /r r grows at low temperatures, whereas the 
opposite is true for x > 0.2, where checkerboard ordering 
is present. 
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FIG. 17: Ratio of diffusion and spin relaxation times as a 
function of temperature at x = 0.02, 0.1, 0.2, 0.3, and 0.4, 
with J/V = 0.2. 

The violation of the DSE relation that we see is much 
smaller in magnitude than is generally seen in experi- 
ments, but this is likely to be because the model we in- 
vestigate displays strong glassy behaviour, whereas ex- 
periments that study the breakdown of the SE and DSE 
relations are usually for fragile glass formers. The exis- 



tence of an extra parameter, the doping x, also allows 
for a richer set of behaviour of T>t rot as a function of 
temperature than is present in supercooled liquids. 



V. DISCUSSION 

We have studied a classical model for a doped antifer- 
romagnet and found signatures of slow and heterogeneous 
dynamics. At low temperatures there is anomalous diffu- 
sion at lengthscales less than l x = l/2yfx and stretched 
exponential spin relaxations. We characterize the tem- 
perature below which the slow spin relaxations occur as 
T ne . For the doping levels we considered T ne is less than 
the higher of the Neel temperature Tjy and the temper- 
ature T c b, below which checkerboard order is observed. 
The timescales associated with both diffusion and spin 
relaxation diverge at low temperatures in an Arrhenius 
fashion, at all doping levels considered, behaviour remi- 
niscent of a strong glass former. Interestingly, we found 
that we were able to fit the temperature dependence over 
the entire temperature range with a sub- Arrhenius form 
for both r r and t s . When we investigated the changes 
in spin configurations between two different times at sev- 
eral different doping levels, we found that at low temper- 
atures, the dynamics were spatially heterogeneous, with 
regions of high and low mobility - we also found that 
spin flips were most likely to occur adjacent to regions 
of high mobility. This spatially heterogeneous dynamics 
also manifests itself in the distributions of local diffusion 
and spin relaxation. The distributions of local square 
displacements indicate that there are fast and slow pop- 
ulations of particles, corresponding respectively, to those 
which have escaped, and those which are trapped in cages 
with size of order l x . The distributions of local spin cor- 
relations also indicate that there are fast and slow sites, 
through the wide range of relaxation times implied by the 
distribution. Finally, we found that when we compare r s 
and r r there are violations of the Debye-Stokes-Einstein 
relation that vary as a function of doping. 

On a first glance at Eq. it might appear strange 
that there are slow dynamics associated with this model. 
There is no quenched disorder, and there is no explicit 
frustration in the interactions as one might expect to see 
in a glassy model. The frustration that leads to slow 
dynamics appears to be hidden in the interplay of anti- 
ferromagnetic interactions and the restrictions on parti- 
cle motion imposed by the level of doping. The clearest 
example of this is the large energy barrier to the mo- 
tion of a hole to the opposite side of a plaquette. One 
way to relate this model to others that have been stud- 
ied previously is to think about integrating out the holes, 
which generates a spin model with plaquette interactions, 
which is a generalized version of the gonihedral models 
that have been studied and found to have suggestions of 
metastable states and glassy dynamics^ 5 - 

Whilst we believe that the study of the model in this 
paper is of interest in itself, we also note that there is 
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a famous class of doped two dimensional antiferromag- 
nets, namely the high T c superconducting cuprates to 
which this work may have some connections. There 
are some very clear differences between our model and 
generic models for these materials, specifically that it is 
classical rather than quantum (and hence cannot show 
superconductivity), and also has Ising spins rather than 
Heisenberg spins, leading to an enhancement of long- 
range order in two dimensions. However, the existence 
of glassy spin and diffusive dynamics, in the absence of 
quenched disorder, along with checkerboard order, even 
in a classical model is very interesting. It also suggests 
that it might be interesting to compare quantities anal- 
ogous to the SE and DSE relations for spin and charge 
relaxations in cuprates. We intend to investigate this 
model with quantum dynamics in future work. This will 
hopefully complement understanding gained in studying 
spin dynamics in doped two dimensional quantum anti- 
ferromagnets with quenched disorder 

The model displays many of the features that one 
would expect in a lattice model for a structural glass (if 
one imagines the spin as corresponding to an orientation 
for a rod-like molecule). The temperature dependence 
of the relaxation times at low temperatures is consistent 
with that of a strong glass-former, although it is inter- 
esting that the temperature dependence over the entire 
temperature range can be fit with a sub-Arrhcnius de- 
pendence. It maybe that dimensionality plays some role 
in the behaviour that we observe - it would be inter- 
esting to study this model in three dimensions to see 



whether similar behaviour is evident there. We note that 
recent experiments on two dimensional glassy films do 
show relaxation times that appear to have sub-Arrhcnius 
behaviour if one fits over a temperature range that strad- 
dles Tg^ 

Another direction of research that we believe may 
be fruitful is, having established the existence of slow 
dynamics in this model to construct a field theoretic 
description, and to compare whether the symmetries 
found in the case of quenched disorder in the Edwards- 
Anderson modeliS are present, and what differences exist 
between the two cases. 
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